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ABSTRACT 


In  order  to  carry  out  the  physics  outlined  in  the  contract 
proposal,  three  major  sets  of  computer  programs  have  been  extensively 
modified  and  extended.  They  calculate  various  electronic  properties  of 
systems  of  atoms  in  the  Unrestricted  Hartree-Fock  (UHF)  approximation. 

The  programs  are:  (1)  a molecular  code  for  finite  clusters  of  atoms,  with 
the  coding  optimized  for  clusters  with  cubic  point  group  symmetry;  (2)  a 
thin  film  code  for  several  interacting  infinite  two-dimensional  sheets  of 
atoms;  and  (3)  a crystal  code  for  an  infinite  number  of  atoms  with  cubic 
symmetry.  These  three  sets  of  programs  provide  a predictive  capability 
for  clusters  and  platelets  of  atoms,  surfaces  and  perfect  crystals 
through  first-principles  calculations  at  the  UHF  level.  The  molecular 
programs  allow  studies  of  imperfect  surfaces  and  crystals  in  which  impurity 
atoms,  dislocations,  vacancies,  surface  reconstructions,  etc. , are  in- 
volved. 

The  molecular  code  is  operational  and  has  been  used  by  us  to 
study  carbon  clusters  containing  five  and  seventeen  atoms  and  fifteen-atom 
lithium  clusters.  Scientists  at  AFAL/DHR  have  used  this  code  to  study  the 
electronic  properties  of  the  nitrogen  impurity  in  diamond.  For  the  carbon 
five-atom  cluster,  it  was  found  that  the  electronic  density  matrix  near  the 
central  atom  is  close  to  that  for  an  infinite  crystal.  The  calculated 
atom-spacing  of  the  seventeen-atom  cluster  is  within  one-percent  of  the 
experimental  value  for  diamond.  This  code  has  also  been  used  in  preliminary 
studies  of  the  electronic  properties  of  small  carbon  platelets. 
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The  current  thin  film  code  uses  s and  p Cartesian  Gaussian  basis 
functions  and  is  designed  for  the  symmetry  of  the  (111)  face  of  diamond. 
The  UHF  code  for  the  perfect  crystal  also  uses  s and  p Cartesian  Gaussian 
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basis  functions. 
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SUMMARY  OF  RESULTS 

| 

I.  Significance  of  the  Research 

Three  distinct  sets  of  Unrestricted  Hartree-Fock  (UHF)  programs 
have  been  developed.  The  molecular  code  (appropriate  for  small  clusters  and 
platelets)  is  fully  operational  (with  an  spd  basis).  The  thin  film  surface 
code  is  operational  (with  an  sp  basis)  for  the  (111)  face  of  diamond.  The 
crystal  program  (with  sp  basis)  has  been  written  but  not  fully  debugged. 

The  thin  film  programs  are  the  only  ones  in  existence  that  have  been  developed 
to  the  UHF  level  of  accuracy.  The  cluster  programs  incorporate  the  fastest 
integral  packages  in  existence,  use  powerful  integral  approximations,  and 
have  been  carefully  optimized  for  the  cubic  point  group  symmetry. 

A large  amount  of  physics  of  direct  interest  to  the  Air  Force  can 


be  done  with  these  programs,  provided  adequate  computer  access  is  available. 

The  strong  similarities  between  the  different  program  sets  (UHF,  Cartesian 
Gaussian  basis  sets,  integral  approximations)  will  allow  valuable  comparisons 
of  the  electronic  properties  of  clusters,  surfaces  and  perfect  crystals. 

These  properties  include  one-electron  energy  levels  and  their  associated  band 
^widths  and  band  gaps,  and  the  atom-projected  density  matrix  (whose  diagonal 
elements  give  the  charge  density  in  the  vicinity  of  any  chosen  atom).  The 
comparison  of  bulk  energy  bands  with  thin  film  energy  bands  is  particularly 
important  for  understanding  the  fundamental  physics  of  various  electronic 
devices.  Cluster  and  platelet  calculations  (which  handle  cases  where  one 
does  not  have  perfect  two  or  three-dimensional  periodicity)  are  very  important 
in  materials  science.  Examples  are  ion  migration  through  crystals,  electron 

, surface  reconstruction,  and  reactions 
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at  surfaces.  Comparisons  of  cluster  calculations  with  bulk  or  film  calcu- 
lations can  test  whether  the  clusters  are  sufficiently  large  to  model  real- 
istically the  physical  phenomena  of  interest.  Thus,  development  of  these 
three  program  sets  is  a necessary  first  step  to  performing  a large  number 
of  accurate  calculations  directly  related  to  practical  problems  in  device 
physics,  materials,  etc. 


II.  Molecular  Code 

The  contract  proposal  states  that  a fast  molecular  code  using 
Cartesian  Gaussians  would  be  employed.  The  two-electron  integral  package 
would  be  extracted  from  the  molecular  code  and  inserted  into  the  surface  and 
bulk  programs.  At  present,  the  fastest  integral  package  for  s and  p Cartesian 
Gaussian  basis  functions  is  contained  in  GAUSSIAN-70  (G-70),  written  at  the 
Carnegie  Institute  of  Technology  by  Professor  Pople  and  his  students.  We 
consequently  chose  G70  to  form  the  nucleus  for  the  present  AFOSR  contract 
work. 

It  was  decided  to  modify  G70  before  extracting  the  relevant  pieces 
for  use  in  the  other  codes.  There  were  several  reasons  for  this  decision. 

(1)  It  was  felt  that  a significant  amount  of  work  relevant  to  the 
contract  could  be  done  with  the  molecular  code  in  cases  where  perfect  two 
or  three  dimensional  translational  symmetry  was  violated.  Examples  are  ion 


migration  in  solids, 


impurity  atom  physics 


in  crystals,  reactions  at  surfaces 


i 


f 


and  atom  reconstruction  at  surfaces. 


(2)  It  was  necessary  to  obtain  considerable  familiarity  with  the 
intricacies  of  the  code  in  order  to  be  able  to  export  sections  to  other 
programs. 

(3)  Most  of  the  programming  algorithms  developed  for  the  molecular 
code  could  be  easily  taken  over  later  to  the  thin  film  and  crystal  codes. 

Consequently,  the  following  major  changes  were  made  to  G70.  (We 
label  the  resulting  code  GH76.) 

(1)  Very  flexible  data  handling  and  program  handling  procedures, 
developed  at  ARL,  were  installed  to  facilitate  program  development. 

(2)  A symmetry  analysis  of  the  molecular  system  reduces  the  number 
of  two-electron  integrals  which  must  be  evaluated  and  handled.  Any  subgroup 
of  the  48  element  cubic  point  group  is  acceptable. 

(3)  Two  integral  approximations  developed  by  us  at  ARL  were 
installed.  One  involves  the  pseudocharge  of  a basis  function  pair,  while 
the  other  involves  the  product  of  two  pseudocharges.  The  first  integral 
approximation  has  proved  extremely  useful  in  cluster  calculations,  while  the 
second  has  saved  so  little  time  that  it  has  been  taken  out  of  GH76.  It  is 
still  extremely  powerful  for  the  surface  and  crystal  programs  however. 

(4)  Routines  which  handle  d basis  functions  were  installed  twice. 
Dr.  P.  C.  Hariharan  had  written  relevant  G70  routines  at  Carnegie  Mellon 
University  for  the  IBM  360.  We  modified  these  routines  for  inclusion  into 
GH76.  Later,  we  heard  that  a new,  extremely  powerful  algorithm  had  been 
developed  for  two-electron  integrals  involving  d basis  functions  by  Professor 
K.  F.  King  and  coworkers  at  the  State  University  of  New  York  at  Buffalo. 
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Their  routines  are  incorporated  into  a molecular  code  called  HONDO.  We 
used  the  appropriate  routines  from  HONDO  (with  their  permission)  in  our 
program  package  (the  H in  GH76) , speeding  up  the  integration  time  for  d 
basis  functions  by  roughly  a factor  of  two. 

(5)  Programs  to  calculate  X-ray  structure  factors  and  directional 
Compton  profiles  were  added  to  GH76.  Because  we  needed  to  study  the  elec- 
tronic density  matrix  in  the  vicinity  of  individual  atoms  (surface  atoms  or 
bulk  atoms),  we  incorporated  the  basic  idea  behind  the  Mulliken  population 
analysis  to  project  out  the  structure  factors,  directional  Compton  profiles 
and  total  energy  for  each  atom  of  the  cluster. 

(6)  Programs  were  written  to  simulate  a cluster  of  atoms  in  an 
external  environment,  i.e.,  to  simulate  the  rest  of  the  crystal  surrounding 
an  impurity  and  its  immediate  neighbors.  The  hope  was  to  calculate  accu- 
rately the  properties  of  a small  cluster  of  atoms,  but  in  an  environment 
which  simulates  the  other  atoms  in  a crystal  or  surface.  We  were  unsuccessful 
in  our  preliminary  attempts,  and  temporarily  abandoned  the  study  in  order  to 
proceed  with  the  development  of  the  thin  film  and  crystal  codes.  However,  we 
feel  that  the  idea  is  sound  and  merits  more  work.  An  alternative  approach 

in  which  many  of  the  electrons  on  peripheral  atoms  are  treated  with  pseudo- 
potentials, while  a more  complete  treatment  is  used  for  the  atoms  in  the 
interior  of  the  cluster,  has  been  examined,  and  shows  considerable  promise. 

(7)  Pseudopotential  methods  were  incorporated  into  our  programs 

so  that  the  atomic  core  electrons  would  not  have  to  be  explicitly  included  in 
the  calculation.  Calculations  are  very  much  faster  when  only  the  valence 
electrons  need  be  considered.  Dr.  Luis  Kahn  of  Battelle  is  a leading  authority 
on  this  procedure.  He  was  kind  enough  to  donate  his  programs,  which  were  then 
appropriately  modified  for  inclusion  in  GH76. 
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(8)  Contour  plotting  and  three-dimensional  perspective  programs 
were  obtained  and  built  into  GH76  so  that  wavefunctions  and  charge  densities 
of  clusters  can  be  graphically  displayed.  This  facility  allows  a visual 
comparison  of  densities  obtained  from  cluster,  surface  and  bulk  calculations. 
It  is  also  a useful  method  for  visually  displaying  the  details  of  the 
electron  bonding  at  surfaces  and  in  crystals. 

(9)  Extensive  modifications  were  made  to  the  portions  of  the 
programs  that  self -consistently  determine  the  UHF  wavefunctions.  The  code 
was  streamlined  and  dynamic  dimensioning  was  introduced  so  that  large 
clusters  could  be  more  easily  run.  A considerably  more  efficient  algorithm 
was  implemented  to  process  the  two-electron  integrals. 

The  following  physics  were  done  with  GH76. 

(1)  Calculations  were  performed  on  a five-atom  carbon  cluster, 
the  atoms  being  arranged  in  a diamond  crystal  structure.  It  was  found  that 
in  the  case  where  all  12  dangling  bond  electrons  have  the  same  spin,  all  of 
the  cluster  results  (band  gap,  valence  band  width,  X-ray  structure  factors, 
directional  Compton  profiles  and  the  binding  energy)  were  close  to  the 
crystalline  results.  Our  Mulliken  projection  procedure  was  necessary  in  this 
investigation  to  project  out  the  density  matrix  for  the  central  atom  of  the 
cluster.  Further  details  are  given  in  Appendix  A. 

(2)  The  total  energy  of  a 17-atom  carbon  cluster  was  calculated 
with  a minimal  basis  set  for  several  atom  spacings.  The  equilibrium  atom 
spacing  was  within  one-percent  of  the  diamond  crystal  spacing.  The  properties 
of  the  central  atom  of  the  cluster  were  again  most  nearly  crystalline  for  the 
case  where  all  36  dangling  bond  electrons  have  the  same  spin. 
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(3)  Fifteen-atom  lithium  crystals  were  calculated  in  collaboration 
with  Dr.  Luis  Kahn.  The  directional  Compton  profiles  for  the  central  atom 
were  close  to  the  crystalline  values.  Further  details  are  given  in  Appendix 
C. 

(4)  A study  was  made  in  collaboration  with  Dr.  Luis  Kahn  to 
determine  the  best  way  to  obtain  the  all-electron  density  matrix  from  the 
pseudopotential  valence-orbital  calculation.  Complications  arise  because 
the  orbitals  obtained  from  the  pseudopotential  calculation  are  not  orthogonal 
to  the  atomic  core  orbitals,  and  do  not  have  the  correct  number  of  nodes. 

We  developed  and  tested  a procedure  which  gives  a density  matrix  within  one- 
percent  of  that  obtained  from  a comparable  (but  much  more  time-consuming) 
all-electron  calculation.  This  knowledge  is  very  important  for  our  surface 
and  crystalline  studies  for  which  the  pseudopotential  procedure  should  be 
used  for  all  of  the  calculations  involving  atoms  heavier  than  neon. 

Further  details  are  given  in  Appendix  B. 

(5)  A study  was  made  comparing  all-electron  and  several  types  of 
pseudopotential  calculations  on  a cluster  of  five  carbon  atoms  in  a diamond 
crystal  structure.  The  study  included  the  use  of  pseudopotentials  for  the 
core  electrons  of  the  peripheral  four  atoms  while  treating  the  central  atom 
fully  (with  all  electrons  included  specifically).  Also,  the  pseudopotential 
calculations  were  carried  out  both  with  normal  (all  electron)  basis  sets  and 
with  valence  basis  sets  developed  specifically  for  pseudopotential  calculations. 
Further  details  are  given  in  Appendix  D. 

The  first  three  results  indicate  that  many  crystalline  properties 
can  be  successfully  obtained  from  small  cluster  calculations  with  G1176  if 
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the  Mulliken  projection  procedure  is  used  to  obtain  the  properties  of  the 
central  atom.  The  fourth  result  shows  that  the  density  matrix  and  associ- 
ated properties  as  well  as  the  energetics  can  be  obtained  to  useful 
precision  with  the  pseudopotential  procedure.  The  fifth  result  demonstrates 
one  practical  procedure  for  including  peripheral  atoms  at  lower  levels  of 
approximation  than  for  interior  atoms  in  cluster  calculations. 

The  GH76  programs  have  been  distributed  in  various  states  of 
completeness  to  several  other  groups.  Dr.  Ronald  Greene  of  AFAL/DHR  used 
them  to  study  the  nitrogen  impurity  in  diamond.  Dr.  J.  Ivey  of  Mound  Lab 
(ERDA)  used  this  program  to  study  the  properties  of  hydrogen  dissolved  in 
metals.  Dr.  P.  Deutch  of  the  Army  Ballistic  Research  Lab  has  used  pieces 
of  GH76,  as  has  Dr.  R.  Jaffe  of  NASA  Ames  Research  Center.  Dr.  T.  C.  Collins 
of  AFOSR  has  used  GH76  to  calculate  the  properties  of  Polymers.  Dr.  R.  J. 
Bartlett  of  Battelle  has  used  these  codes  for  a study  of  the  electronic 
structure  of  catalytic  nickel  surfaces.  The  program  has  also  been  provided 
to  Dr.  M.  W.  Ribarsky  of  the  Georgia  Institute  of  Technology  and  to  Dr.  D.  W. 
Jennison  of  Sandia  Laboratories. 

III.  Thin  Film  and  Crystal  Codes 

It  became  a major  programming  effort  to  incorporate  the  GH76  s-p 
integral  package  into  the  surface  and  bulk  programs.  Previously,  both  codes 
handled  each  two-electron  integral  over  four  basis  functions  separately. 

GH76  gains  its  speed  and  efficiency  from  doing  blocks  of  integrals  simutaneous ly. 
P-orbitals  have  the  same  Gaussian  exponents  as  s orbitals,  and  all  of  the 
s-x-y-z  blocks  of  integrals  for  all  four  functions  are  calculated  at  once. 
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The  incorporation  of  this  powerful,  time-saving  feature  involved  substantial 
changes  in  the  logical  structure  of  the  programs.  Consequently,  what 
initially  appeared  at  contract-proposal  time  to  require  some  weeks  of  work 
turned  out  to  require  many  man-months  of  work.  Long  delays  and  limitations 
on  access  to  the  Wright-Patterson  Air  Force  Base  computer  hampered  this  work 
considerably  and  limited  its  scope. 

The  thin  film  programs  were  converted  first.  They  were  debugged 
as  much  as  possible  on  the  Battelle  computers  at  Battelle  expense.  The  pro- 
grams were  then  run  on  the  WPAFB  computers  with  a realistic  Gaussian  basis. 
Early  results  appeared  very  reasonable  for  a two-atom  layer  film  of  diamond, 
but  extensive  production  computations  were  impossible  due  to  very  limited 
computer  access  and  very  slow  service  at  WPAFB. 

The  crystal  programs  have  also  been  restructured  and  rewritten  to 
handle  integrals  in  blocks.  They  have  first  been  debugged  for  small  calcu- 
lations on  the  Battelle  computer  at  Battelle  expense  (with  a trial  basis), 
but  final  debugging  with  a more  realistic  basis  at  WPAFB  could  not  be  com- 
pleted due  to  computer  access  limitations. 

Changes  to  the  surface  and  crystal  programs  included  the  following: 

(1)  The  GH76  one-  and  two-electron  integral  packages  for  s-p  func- 
tions were  installed,  with  the  appropriate  logical  restructuring. 

(2)  The  ARL  Gaussian  lobe-oriented  integral  approximation  was 
scrapped  and  replaced  by  an  approximation  based  upon  a moment  expansion  of 
the  charge  distribution  of  the  basis  function  product. 

(3)  The  Ewald  procedure  (which  finds  the  potential  of  the  infi- 
nitely repeating  array  of  charge  distributions)  was  revised  to  use  this  same 
moment  expansion. 

L _____ 
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Extensions  to  include  d basis  functions  were  contemplated,  but 
implementation  could  not  be  carried  out  due  to  time  limitations  and  computer 
access  limitations. 
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Abstract.  X-ray  structure  factors  and  directional  Compton  profiles  are  calculated  for  diamond 
(cubic  carbon)  using  a self-consistent  field  unrestricted  Hartree-Fock  procedure  on  a cluster 
of  five  carbon  atoms.  A Mulliken  population  analysis  procedure  is  used  to  isolate  results  to 
be  associated  with  the  central  atom.  The  choice  of  a high  spin  state  for  the  cluster  yields 
reasonable  agreement  with  the  experimental  and  theoretical  results  for  crystalline  diamond. 


1.  Introduction 

. 

The  Hartree-Fock  (HF)  equations  have  been  solved  previously  for  diamond , and 
selected  crystalline  properties  have  been  compared  with  experiment.  The  calculated 
x-ray  structure  factors  (Euwema  et  al  1973)  and  directional  Compton  profiles  (Wepfer 
' et  al  1 974)  which  r leasure  the  electronic  charge  and  momentum  distributions  respectively, 

are  in  satisfactory  agreement  with  experimental  results  (Clark  1964,  quoting  Gotlicher 
and  Wolfel  1959;  Reed  and  Eisenberger  1972).  The  one-electron  energies  and  the  co- 
hesive energy  are  what  one  would  expect  when  correlation  has  been  completely  neglected. 
However,  the  crystalline  HF  formalism  has  the  disadvantages  of  requiring  much  computer 
time  and  effort,  and  of  requiring  complete  translational  symmetry.  This  formalism  and 
the  present  crystalline  HF  computer  programs  cannot  be  readily  extended  to  treat 
more  complex  crystals,  impurities,  or  surfaces. 

The  purpose  of  this  paper  is  to  demonstrate  the  feasibility  of  obtaining  certain  crystall- 
ine results  from  unrestricted  HF  (UHF)  (Pople  and  Nesbet  1954)  calculations  on  very 
small  clusters.  In  particular,  we  shall  display  x-ray  structure  factors  and  directional 
Compton  profiles  for  the  central  atom  of  a five-atom  diamond  cluster  as  well  as  the 
average  over  the  cluster.  Results  for  the  central  atom  agree  reasonably  well  with  those 
of  the  infinite  crystal.  In  order  to  obtain  this  agreement,  it  was  necessary  to  choose  the 
correct  spin  multiplicity  for  the  cluster,  and  to  employ  the  Mulliken  population  analysis 
procedure  to  divide  the  density  matrix  into  portions  associated  with  each  of  the  indivi- 
dual atoms.  The  spin  multiplicity  results  have  important  implications  for  surface  calcula- 
tions on  covalently  bonded  materials. 

In  §2,  wc  discuss  the  computational  details,  including  the  Mulliken  prescription  used 
to  divide  up  the  density  matrix.  The  question  of  spin  multiplicity  is  treated  in  §3,  while 
§4  and  §5  are  devoted  to  discussions  of  the  structure  factors  and  Compton  profiles.  The 
results  are  then  summarized  in  §6. 

fPrescnt  address:  Chemical  Engineering  Division,  Argonnc  National  Laboratory,  Argonnc,  Illinois 60439,  USA. 
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2.  Computational  method 

The  self-consistent-field  UHF  procedure  has  been  used  to  calculate  one-electron  energies, 
total  energies,  and  one-electron  orbitals,  for  a five-atom  carbon  cluster,  employing 
a modified  version  of  the  Gaussian  70  program  package  of  Pople  el  al. t We  used  the 
1 ls/5p  cartesian  gaussian  carbon  basis  of  Pople  and  Binkley  1 1975).  contracted  to  4s/3p. 
The  carbon  atoms  have  the  coordinates  (0,  0,  0),  a/M  I.  1,  l),a/'4(l,  —1.  — 1 ),  ei/4<  — 1 . —I. 
1)  and  a/4(-l.  1,  --1).  The  diamond  lattice  constant,  a,  is  3-56  A.  The  atoms  form  a 
tetrahedron,  and  the  distance  between  the  centre  atom  and  any  one  of  those  at  the  corners 
is  the  nearest-neighbour  distance  in  diamond. 

The  first-order  density  matrix  is  obtained  by  summing  over  the  filled  LIHF  orbitals 

p(r,r')  = V 4'M(r)'P’(r').  (1) 

filled 

state* 

The  Mulliken  population  analysis  procedure  (Mulliken  1955)  is  then  used  to  assign  parts 
of  p to  the  individual  atoms  in  the  cluster  as  follows.  The  orbitals  are  expanded  in  the 
contracted  gaussian  basis  functions,  <)>Jr  — Ra),  associated  with  the  atoms,  {a at 
positions  {RJ.  Expanding  the  density  matrix  in  terms  of  the  basis  functions,  one  obtains 

f Hr,  r')  = T £ Kl  <t>Jr  - Ra)  4>„(r  - RJ  (2 a) 

ah  afi 

= £ pjr  »•')  (2b> 

ab 

where  pab  is  that  part  of  the  density  matrix  associated  with  the  pair  of  atoms  a and  h. 
Following  the  Mulliken  prescription,  we  assign  half  of  each  term.  pjl(.  to  atomic  site  a 
and  half  to  site  h : 

p(r,  r')  = £ [pcc  (r,  r')  + -|  V (pc<f  (r,  r')  + pjc(r,  r'))]  (3a) 

c d*c 

= Tpf(r,r').  (3b) 

C 

This  decomposition  of  p can  then  be  used  to  calculate  the  structure  factors  and  Compton 
profiles  for  the  individual  atoms  of  the  cluster. 

The  Mulliken  decomposition  is  basis  dependent  for  small  molecules.  An  extreme 
example  would  involve  using  a rich  basis  on  the  central  atom  of  the  five-atom  cluster  and 
no  basis  functions  on  the  outer  four  atoms.  In  this  case,  the  central  atom  would  be  assigned 
the  entire  density  matrix.  However,  in  an  infinite  crystal,  any  appropriate  basis  would 
presumably  have  the  full  translational  and  rotational  symmetry  of  the  crystal's  space 
group  In  this  limit,  every  equivalent  unit  of  the  crystal  is  assigned  an  equal  portion  of 
the  density  matrix.  Thus,  to  the  extent  that  the  unit  at  the  centre  of  the  cluster  is  far  from 
the  boundaries,  the  Mulliken  decomposition  should  give  crystalline-type  .esults  lor  the 
central  unit. 


3.  Spin  multiplicity 

The  correct  many-body  wavcfunction  for  the  carbon  cluster  may  be  written  as  a sum 

t The  Gaussian  70  proprafn  package  was  developed  at  the  Carnegie  Mellon  University  by  Professor  J A Pople 
and  his  students  (W  J llcdrc.  W A Lathan.  R Ditchficld.  M I)  Newton  and  J A Pople,  Quantum  Chemistry 
Program  Exchange.  Indiana  University,  Program  236). 
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over  Slater  determinants.  In  the  UHF  approximation,  however,  one  selects  the  single 
Slater  determinant  that  best  represents  the  system  under  consideration.  In  particular,  one 
must  specify  the  spint.  S.  or  equivalently,  the  multiplicity.  2 S'  4-  I.  of  the  Slater  determi- 
nant. The  spin  determines  the  occupancy  of  the  spin-up  and  spin-down  orbitals.  For  the 
five-atom  carbon  cluster,  (15  + S)  spin-up  states  and  (15  - S)  spin-down  states  are 
occupied. 

To  obtain  realistic  surface  states,  one  wants  the  spin  state  that  gives  the  lowest  energy 
for  the  cluster.  Nothing  guarantees,  however,  that  the  lowest  energy  state  will  give  the 
best  crystalline  properties  at  the  central  atom.  To  study  these  equations,  UHF  calcula- 
tions have  been  performed  for  cluster  spins  ranging  from  0 to  6.  Selected  results  are  shown 
in  table  1 and  in  figure  1. 


Tablet.  Dependence  upon  cluster  spin  (S|  of  the  total  cluster  energy  (te).  charge  on  the 
central  atom  (Q).  valence  band  width  (vbw),  band  gap  (H(.)  and  the  difference  between  the 
(100)  and  (111)  directional  Compton  profiles  at  q = 0(A7(0)).  The  row  labeled  'atom'  gives 
the  total  energy  of  five  non-interacting  carbon  atoms  with  the  same  gaussian  basis  The  row 
labeled  ‘crystal'  gives  HF  crystal  results.  The  labels  ILT)  and  (gt)  indicate  whether  the 
cluster  results  are  expected  to  be  ‘less  than'  or  greater  than'  the  crystalline  or  atomic  results. 
AU  energies  are  in  au. 


s 

TE 

Q 

VBW 

BCi 

A7(0) 

atoms 

(IT) 

- 188-4230 

6 

crystal 

KIT) 

- 189  3348 

6 

(ET)  Ml 

(gt)  0-52 

0 10 

0 

- 1884172 

6 70 

115 

009 

-0-22 

1 

- 188  3169 

6-75 

1-20 

0-10 

-0-20 

2 

-188  4668 

6 85 

1 16 

033 

-019  ' 

3 

- 188  3789 

662 

1-14 

035 

-0  12 

4 

- 188-2701 

6-31 

1-03 

046 

-0-02 

5 

-188-4153 

6-03 

091 

0-57 

0-10 

6 

-188  6392 

6 17 

0-78 

067 

0-05 

The  total  cluster  energy  dips  at  spin  2,  but  is  lowest  for  spin  6 (by  6 cV  over  the  spin  0 
energy).  From  the  figure  we  see  that  only  the  clusters  with  spins  2 and  6 are  bound  with 
respect  to  the  energy  of  the  free  atoms.  The  charge  on  the  central  atom  will  be  six  electrons 
in  an  infinite  crystal.  The  central  atom's  charge  for  the  cluster  is  closest  to  this  value  of 
six  for  spins  5 and  6.  For  the  central  atom,  the  difference  between  the  directional  Compton 
profiles  in  the  (100)  and  (111)  directions  has  the  wrong  sign  at  low  spin,  and  is  closest  to 
the  crystalline  result  for  spins  5 and  6.  .lecause  the  one-electron  eigenvalues  for  the 
cluster  form  a subset  of  the  crystalline  eigenvalues,  the  valence  band  width  (energy  spread 
of  the  non-ls-likc  filled  orbitals)  of  the  cluster  should  be  narrower  than  the  crystalline 
HF  bandwidth.  Similarly,  the  bandgap  (energy  difference  between  the  highest  occupied 
and  lowest  empty  orbitals)  between  empty  and  occupied  states  should  be  larger  for  the 
cluster  than  for  the  crystal.  These  inequalities  also  are  met  only  at  high  spin. 

We  conclude  for  these  reasons  that  the  spin  6 state  of  the  cluster  gives  the  lowest 
total  energy,  the  most  realistic  surface  states,  and  the  most  crystalline  results  for  the 

t In  UIIF  calculations,  the  wavefunction  is  not  necessarily  an  eigenfunction  of  the  spin  However,  in  the  case 
of  C5,  the  expectation  values  of  the  spin  were  close  to  the  integer  values  reported  here 
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central  atom.  In  this  state,  21  electrons  have  spin  up  and  Q have  spin  down.  Five  electrons 
of  each  spin  are  in  carbon  Is  core  states.  Four  of  each  spin  are  in  bonding  orbitals  between 
the  central  atom  and  its  nearest  neighbours.  The  remaining  12  electrons  all  have  spin 


Figure!.  The  total  energy  of  the  five-atom  cluster  carbon  as  a function  of  the  spin.  S.  (more 
accurately.  At,).  The  broken  line  gives  the  energy  of  five  isolated  carbon  atoms.  The  value  for 
the  infinite  HF  crystal  would  be  - 189-3348  au. 

up.  They  are  therefore  spatially  orthogonal  to  each  other  and  to  the  central  bonding 
orbitals,  so  that  they  point  roughly  towards  the  (missing)  12  carbon  atoms  at  the  next- 
nearest-neighbour  positions  in  diamond:  a/2(  + l,  +1,  0),  a/2(+l,  0+1)  and 
a/ 2(0,  ± 1,  ± 1).  They  have  gained  6eV  of  exchange  energy  over  the  spin  0 state  of  the 
cluster. 

This  is  an  example  of  Hund’s  rule  (Herzberg  1944),  extended  to  molecules:  ‘Of  the 
configurations  given  by  equivalent  electrons,  those  with  the  greatest  multiplicity  lie 
deepest . . If  this  rule  were  to  hold  for  a more  accurate  representation  of  a surface  on  a 
covalently  bound  material,  the  electrons  in  the  dangling  bond  surface  states  would  all 
have  the  same  spin. 


4.  X-ray  structure  factors 

X-ray  structure  factors  give  a quantitative  measure  of  the  electronic  charge  distribution 
by  providing  the  Fourier  transform  of  the  charge  density  for  a given  reciprocal  lattice 
vector  of  the  crystal. 

S(k)  = fcxp(ikr)p(r,r)d(r).  (4) 

In  order  to  compare  the  structure  factors  for  each  atom  of  the  cluster  with  those  of 
the  infinite  crystal,  we  must  first  multiply  the  atom’s  structure  factor  obtained  from  the 
Mullikcn-decomposcd  density  matrix  by  exp {--k-Rj  to  remove  a meaningless  phase, 
and  then  impose  the  inversion  symmetry  of  the  diamond  cubic  crystal  structure. 
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P(r , r ')  = p(T~  r,T-  r ') 
for 

T = a/4  (1, 1,  1). 


(5a) 

(Sb) 


We  thereby  form  for  an  atom,  a,  of  the  cluster 

Sjk)  = S„(M  + exp  (ik-T)Sf  — k).  (6) 

The  results  are  shown  in  table  2.  With  the  exception  of  the  (222)  direction,  the  crystalline 
HF  and  the  experimental  structure  factors  are  within  2%  of  each  other,  which  is  within 
both  experimental  and  computational  uncertainty.  The  tendency  of  the  (222)  results  to 
be  too  small  has  been  seen  in  other  materials  as  well  as  diamond,  and  is  probably  due  to 
a lack  of  d functions  in  the  carbon  basis.  Results  for  the  central  atom  of  the  cluster  are 
considerably  closer  to  the  crystalline  results  than  are  the  results  for  the  surface  atoms  or 
for  the  average  over  all  cluster  atoms.  Even  for  a cluster  of  as  few  as  five  atoms,  the  charge 
distribution  associated  with  the  central  atom  is  already  close  to  that  of  an  atom  in  the 
perfect  crystal. 


Table  2.  X-ray  structure  factors  from  experiment,  from  crystal  HF,  for  the  five  individual 
carbon  atoms  at  the  given  coordinates,  for  the  average  over  the  surface  atoms,  and  for  the 
average  over  the  cluster.  Diamond  inversion  symmetry  is  imposed  as  is  described  in  the  text. 


000 

111 

220 

311 

222 

Experiment 

6 

332 

1-98 

166 

0-14 

Crystal  HF 

6 

3-27 

1-94 

1-69 

0-08 

(0,  0,  0) 

6-16 

3-21 

192 

1-68 

0-07 

a/40,  1,  1) 

596 

322 

1-99 

1 85 

001 

o/4(l,  -1,-1) 

5-96 

2-99 

1-97 

1-74 

002 

fl/4(  — 1.  1,  -1) 

596 

299 

1-97 

1-76 

002 

a/4(—  1,  -1,  1) 

596 

299 

1 99 

1-76 

0-02 

Surface 

596 

304 

1-98 

1-78 

0-02 

Cluster 

6 

308 

197 

1-76 

003 

5.  Directional  Compton  profiles 

Compton  profiles  are  related  to  the  electron  momentum  distribution,  denoted  t(k).  In 
the  HF  approximation, 

t(k)  = £ |«jH*)|2  =,Ar  fcxP  Cifc '(r  ~ r')]  r') d(r) d(r').  (7) 

filled  J 

stales 

In  particular,  in  the  impulse  approximation,  the  Compton  profile  in  the  direction  n 
is  given  by 

Jm{  </)  = J t(k)JUi  - k n)  d(k).  (8) 

This  is  just  an  integral  of  the  momentum  distribution  over  a plane  in  k space  at  a distance. 
q.  from  the  origin  taken  in  the  direction  of  n.  The  dependence  of  J.{q)  on  the  direction,  n. 
thus  gives  information  about  the  bonding  orbitals  of  diamond.  Additionally,  each 
Compton  profile  satisfies  the  sum  rule 

j JJq)  dq  = \ t (k)  d(lc)  = number  of  electrons. 


(9) 
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It  is  well  known  that  in  materials  like  diamond,  the  crystal  wavefunction  has  more 
high-momentum  (high-k)  components  that  it  does  in  the  atom,  since  the  crystalline 
electrons  are  localized  into  a smaller  volume.  Therefore,  because  of  the  sum  rule,  the 
Compton  profiles  at  zero  momentum  (q  = 0)  fall  in  going  from  the  isolated  atom  to  the 
crystal.  This  is  illustrated  in  figure  2,  where  we  see  the  progressive  decrease  in  J, , ,(q) 


Figure  1 Compton  profiles  in  the  ( 1 1 1 ) direction  for  the  carbon  atomf ),  Cj  cluster 

( ),  C5  cluster  ( ) and  HF  crystal  I ). 

from  the  atom  to  the  C2  molecule,  to  the  C5  cluster,  and  then  to  the  infinite  crystal. 
Because  the  cluster  average  for  all  five  atoms  is  very  close  to  the  value  for  the  central 
atom,  we  show  the  Compton  profile  for  the  average  in  the  figure.  It  is  obvious  from  the 
figure  that  a substantial  amount  of  high-momentum  character  is  still  missing  from  the 
cluster  Compton  profile. 

Figures  3(a)  (</)  show  the  differences  between  the  Compton  profiles  in  various  direc- 
tions. These  differences  are  in  satisfactory  agreement  with  experiment  and  with  the  crystal 
HF  calculations  for  all  of  the  cases  except  the  (100)  - (110).  Significantly,  the  (110) 
direction  in  diamond  corresponds  to  the  direction  of  the  next-nearest-neighbour  of  the 
central  atom,  and  these  atoms  are  not  included  in  the  cluster. 


6.  Summary  and  conclusions 

We  have  seen  that  it  is  possible  to  obtain  fairly  reasonable  crystalline  x-ray  structure 
factors  and  directional  Compton  profile  differences  for  diamond  from  a cluster  containing 
as  few  as  five  atoms.  This  cluster  includes  the  dominant  interactions  in  such  a solid  of 
the  central  atom  with  its  nearest  neighbours.  It  is  essential,  however,  to  choose  that  spin 
multiplicity  for  which  all  of  the  dangling  orbitals  have  the  same  spin.  One  can  then 
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eliminate  the  effects  of  these  surface  states  from  results  to  be  attributed  to  the  bulk  by 
using  the  Mulliken  population  analysis  which  partitions  the  density  matrix  with  respect 
to  the  individual  atoms.  From  the  results  given  above  for  five  atoms,  it  is  suggested 
that  in  the  case  of  a cluster  of  17  atoms  (one  carbon  atom  plus  both  its  nearest  and  next- 
nearest  neighbours)  a spin  state  of  18  (36  dangling  bonds  at  the  surface)  might  give  results 
for  the  central  atom  which  would  be  in  quite  satisfactory  agreement  with  crystalline 


FiRurcJ.  Directional  Compton  profile  differences  for  experiment  (•  • •(.  crystalline  ItF 

I X cluster  average  ( X and  central  atom  of  the  cluster  ( X for  (u) 

(100)  - (HlX(fc)flOO  - (llOX(c)flOO)  - (I  I2X  M)  (100)  - (221). 


results.  Further  work  is  also  in  progress  to  investigate  the  relationship  between  the  critical 
dependence  of  the  cluster  results  and  the  spin  multiplicity,  and  more  realistic  models 
of  the  diamond  surface. 
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The  accuracy  of  first-order  density  matrices  obtained  by 
electronic  structure  methods  involving  the  valence  electrons 
only 
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First-order  density  matrices  are  obtained  from  all-electron  unrestricted-Hartree-Fock  (UHF)  calculations 
and  from  elTective-core-potential  valence-electron  UHF  calculations  for  the  silicon  and  argon  atoms,  and 
for  the  Si,  molecule.  The  density  matrices  are  compared  by  means  of  their  derived  X-ray  scattering  factors, 
directional  Compton  profiles,  and  the  charge  density,  at  selected  points  in  k and  r space. 


I.  INTRODUCTION 

A central  problem  in  molecular  and  solid-state  com- 
putational physics  involves  the  calculation  from  first 
principles  of  the  electronic  first-order  density  matrix.1'2 
Quantities  derivable  from  the  density  matrix,  such 
as  X-ray  scattering  factors  and  directional  Compton  pro- 
files are  measurable  to  high  precision— often  to  within 
one  percent. 3 

For  molecular  systems  with  constituent  atoms  heavier 
than  neon,  it  rapidly  becomes  very  expensive  to  perform 
electronic  structure  calculations  involving  all  of  the 
electrons.  Procedures  based  upon  pseudopotential  the- 
ory have  recently  been  developed  for  reducing  the  all- 
electron (AE)  molecular  problem  to  a valence -electron 
(VE)  problem. 4 The  tightly  bound  core  electrons 
(1s22sz2/>6  for  silicon  and  argon,  for  example)  are  as- 
sumed to  have  the  same  wavefunctions  as  in  the  isolated 
atoms.  An  effective  core  potential  (ECP)  is  constructed 
directly  from  the  atomic  Hartree-Fock  equations,  fol- 
lowing a well  defined  pseudo-orbital  transformation,  to 
represent  the  effect  of  the  frozen  core  orbitals  on  the 
remaining  electrons.  The  VE  problem  in  the  field  of 
the  effective  core  potentials  is  then  solved  by  the  stan- 
dard electronic  structure  methods. 

The  purpose  of  this  paper  is  to  compare  first-order 
density  matrices  calculated  from  AE  UHF  calculations 
with  those  obtained  from  ECP-VE  UHF  calculations  to 
obtain  an  estimate  of  the  percent  deviation,  and  hence 
of  the  percent  reliability  of  the  ECP  density  matrix. 
Comparisons  will  be  made  of  three  quantities  derived 
from  the  density  matrix:  X-ray  scattering  factors, 
directional  Compton  profiles,  and  the  electronic  charge 
density. 

We  have  chosen  three  systems  for  comparison:  the 
open-shell  silicon  atom,  the  closed  shell  argon  atom, 
and  the  diatomic  molecule  Si2. 


and  the  expectation  value  of  any  one-electron  operator 
is  given  by 

(P>=  f 6(r'  - r)P(r')p(r,r')c/(r)(/(r')  . (2) 

The  all- electron  wavefunction  is 

r„) 

= A[<pu  a--  - <P„aa<PuP---  <pnbf30v(l, , N-  2n)] , (3) 

where  {<fiia<plb,  i=  1, ...,«}  are  the  core  orbitals,  and 
0„  is  the  valence  wavefunction.  In  the  UHF  approxima- 
tion, the  valence  wavefunction  is  given  by 

©i>=  <P(n*l)aa  • • * V (n*na)aa(P (n.l)bP  • ' * (rent.  160  (4) 

(2 n + na  + nb=N)  . 

If  the  valence  and  core  orbitals  form  an  orthogonal  set, 
then  the  first-order  density  matrix  is  given  by  the  well- 
known  expression 

n»na 

p(r,  r')=  <Pk.b)(pfa(r')+  H <pti( r)<pt„(r') . (5) 

*=i  *-i 

However,  one  of  the  essential  changes  introduced  in  the 
reduction  of  the  AE  problem  to  a VE  problem  is  the  lift- 
ing of  the  core  orthogonality  requirement  on  the  valence 
orbitals.4  The  orthogonality  effects  are  folded  into  the 
effective  core  potentials.  The  modified  valence  orbit- 
als are  smooth  nodeless  orbitals4  which  are  not  orthog- 
onal to  the  core  orbitals.  It  is  therefore  incorrect  to 
use  Eq.  (5)  for  the  density  matrix,  when  the  valence  or- 
bitals are  determined  from  a VE  calculation. 

In  order  to  obtain  the  correct  expression  for  the  one- 
electron  density  matrix,  it  is  convenient  to  transform 
from  the  nonorthogonal  sets  of  core  and  valence  orbitals, 
{<Pka>  + na\  and  {<pbb,  k = 1, . . . , n + nb}  to  sets 

of  orthonormal  orbitals,  {i pkt,  k = l ,...,»  + >w}  and  {(//,>, 

k=  1 H + wfij.  The  orthonormalizing  transformation 

may  generally  be  written  as 

n+na 


II.  CALCULATION  OF  THE  ONE-ELECTRON  DENSITY 
MATRIX 

4>k.=  kL  <PlaAt>  k=l,n+na  , 

<•1 

H+flfe 

(6a) 

The  one-electron  density  matrix  is  obtained  from  the 
electronic  wavefunction  by  the  relation5 

</'*»=  53  <Pu *=1  ,»*nb  , 

<-i 

where  the  A and  D matrices  satisfy  the  conditions 

(6b) 

p(t,v')=N  J[*(r,rz r„) 

A-AaA  = 1 , 

(7a) 

x +*(r',  r2, . . . , r„)d(r2). . .r/(r/v)/('H'>],  (1) 

B‘A„B  1 , 

(7b) 
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and  where 

and 


~ (Vial  Vjt) 

We  note  that  the  UHF  wavefunction  constructed  from 
the  nonorthogonal  orbitals  is,  apart  from  a constant, 
equal  to  a UHF  wavefunction  constructed  from  the 
orthogonal  orbitals.  Therefore,  the  correct  expres- 
sion for  the  one-electron  density  matrix  is 

p(r , r')=  £ £ ^.(r)^(r')  . (8) 

*■1  *«1 

Inserting  Eq.  (6)  in  Eq.  (8),  we  obtain 

n*na  »»»na 


p(r,r')=  £ £v>i.(r)(AA‘)(i^<i(r') 

l-l  i-i 
n*nb  n*nb 

+ £ £*„(r)(BB%vjl(r'). 


ci  j- 1 


(9) 


The  general  form  of  the  A and  B transformation  ma- 
trices satisfying  Eqs.  (7)  is* 


a = a;i/2ua  , 
B=A-.1/2Ub  , 


(10a) 

(10b) 


where  UA  and  UB  are  arbitrary  unitary  matrices.  It 
follows  readily  that 


aa»  = a;‘ 


(1  la) 


and 


BB*=  a; 


(lib) 


Therefore,  the  correct  expression  for  the  UHF  one- 
electron  density  matrix  is 


rr+nq  n+na 


P(r,r’)=  £ £ ^JrXA;1),, <*>/>') 

<■1  }•  1 

n»n6  n*nb 

+ £ £ vJr)U-Bl)u<p%(r')  . 


(12) 


TABLE  I.  All-electron  (AE)  and  ECP  valence  electron  eigen- 
values for  argon  and  silicon  atoms  and  for  the  molecule  SU. 

The  double  entries  for  atomic  silicon  give  the  spin-up  and  spin- 
down  llllF  Eigenvalues.* 


AE 

ECP-VE 

Argon 

3s 

-1.275 

-1.274 

3 P 

-0.589 

-0.589 

Silicon 

3s 

-0.618  -0.461 

-0.019  -0.458 

3 P,., 

-0.297 

-0.297 

Si2 

4(Tj 

-0.690 

-0.683 

4a, 

-0.489 

-0.493 

-0.281 

-0.277 

•All  quantities  are  expressed  in  atomic  units.  See  Ref.  9. 


(ii)  the  X-ray  scattering  factor  S(k)  = /a(r)e“  rd(r),  which 
Is  the  Fourier  transform  of  the  charge  density,  and  (iii) 
the  directional  Compton  profile  (in  the  impulse  approxi- 
mation) 


Jn(q)=  j f(k)6(n- k-q)d(k)  , 


To  obtain  the  above,  we  have  used  some  special  proper- 
ties of  the  UHF  wavefunction.  For  multideterminantal 
wavefunctions,  <i>,  the  transformations  of  the  type  (6) 
will,  in  general,  transform  into  a manifold  of  deter- 
minants of  types  possibly  not  present  previously.  How- 
ever, once  is  reexpressed  in  terms  of  the  orthonormal 
set  of  orbitals,  it  is  straightforward  to  obtain  the  one- 
electron  density  matrix  by  use  of  Eq.  (1). 

The  procedure  outlined  above  for  constructing  the 
density  matrix  from  nonorthogonal  core  and  valence  or- 
bitals [Eq.  (12))  can  be  contrasted  with  a method  out- 
lined earlier,  Eq.  (9)  of  Ref.  7,  in  which  the  valence 
orbitals  were  explicitly  orthogonalized  to  the  core  or- 
bitals. However,  in  Ref.  7,  the  resulting  core-or- 
thogonalized  valence  orbitals  were  never  made  mutually 
orthogonal  to  each  other. 

III.  ONE-ELECTRON  PROPERTIES  CALCULATED 
FROM  THE  DENSITY  MATRIX 

The  one-electron  properties  considered  here  are  (1) 
the  charge  density  o(r)  =/Mr’  - r)p(r,  r)rf(r')  = p(r,  r), 
which  Is  simply  the  diagonal  part  of  the  density  matrix, 


which  is  obtained  by  integrating  the  electron  momentum 
distribution,  f(k),  over  a plane  normal  to  a vector  n at  a 
distance  q from  the  origin  in  momentum  space.  The 
momentum  distribution,  /( k),  is  obtained  by  separate 
Fourier  transforms  in  r and  r' of  p(r,  r'), 

<(k)=  (2xF  / /«“'rp(r,r')e-«  r'd(r)d(r')  . 

IV.  CALCULATIONAL  DETAILS 

The  first  system  chosen  for  detailed  study  is  the  open 
shell  silicon  atom  in  the  3P  state.  The  valence  spin- 
orbital  product  is 

(3s„ ot)(3stP)(3p„ o,)(3/)yo  a)  . 

The  basis  consists  of  Huzinaga’s  (14sl0/>)  Cartesian 
Gaussian  basis  contracted  to  a [4s3/>]  set.®  The  total 
energy  for  this  basis  is  - 288.  8465  hartrees. 9 

The  closed  shell  argon  atom  was  chosen  as  the  second 
system.  We  use  Huzinaga’s  (l4sl0/>)  Gaussian  basis  for 
argon,®  contracted  to  a [4s3/>]  set.  The  total  energy  is 
then  - 526.  79957  hartrees. 

The  same  silicon  Gaussian  basis  was  used  for  the 
UHF  calculation  on  (he  lZ*  state  of  Si2.  The  interatomic 
separation  used  is  4.4352  bohr.  9 The  total  energy  is 
then  - 577.6263  hartrees.  The  valence  configuration  is 

(4oJ)(4o2)(27r^,)(2  wfB)  . 

In  all  cases,  the  same  basis  was  used  for  the  AE  cal- 
culation and  for  the  ECP-VE  calculation,  although  in 
practice  a much  smaller  basis  suffices  for  the  ECP-VE 
calculation  because  of  the  smoothness  of  the  valence 
pseudo- orbitals. 4 Table  I shows  a conqiarison  of  the 
orbital  energies  from  the  AE  and  ECP-VE  calculations. 
This  comparison  shows  agreement  of  the  eigenvalues  to 
within  0.007  hartrees. 

The  parameters  specifying  the  effective  core  poten- 
tials  for  both  silicon  and  argon  were  taken  from  the  re- 
cent work  by  Topiol. 10 
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TABLE  II.  Calculated  X-ray  structure  factors,  S(/C),  and 
directional  compton  profiles,  </„(</),  for  atomic  silicon. *,b 


k or  q 

AE 

Eq.  (12)(%)c 

Eq.  (5)  (%)" 

SIOO 

0.2 

13.77 

0.00 

0.06 

0.4 

13. 14 

0.01 

0.22 

1 

10.55 

0.03 

0.  84 

2 

8.11 

0.13 

0.29 

4 

5.29 

0.  38 

-1.  74 

S0oi 

0.2 

13.  84 

0.00 

0.04 

0.4 

13.38 

0.01 

0.  16 

1 

11.20 

0.03 

0.69 

2 

8.42 

0.  13 

0.59 

4 

5.  28 

0.39 

-1.38 

0 

4.74 

-0.09 

2. 19 

0.4 

3.  86 

-0.  10 

2.41 

1 

1.73 

-0. 18 

1.05 

2 

1.04 

-0. 12 

-3.  53 

4 

0.37 

0.62 

-1.97 

•Tool 

0 

5.83 

-0.06 

1.80 

0.4 

3.70 

-0.09 

3.05 

1 

1.58 

-0.23 

0.11 

2 

1.05 

-0.14 

-4.00 

4 

0.37 

0.63 

-1.25 

“The  3P  silicon  wavefunction  has  the  3p,  and  3py  one-electron 
orbitals  occupied.  Thus  the  SIOO  and  </100  refer  to  K and  n 
within  the  plane  of  symmetry. 

bAll  quantities  are  expressed  in  atomic  units.  See  Ref.  9. 

'The  column  headed  “Eq.  12"  gives  the  percent  deviation  of 
ECP  results  from  the  all-electron  results  (AE),  when  Eq.  (12) 
Is  used  to  calculate  the  total  (core  plus  valence)  density 
matrix. 

"The  column  headed  “Eq.  (5)"  is  similarly  defined. 


V.  RESULTS 

The  all-electron  results  and  the  precent  deviations  in 
the  ECP  valence  electron  results,  by  use  of  both  the 
approximate  density  matrix  expression,  Eq.  (5),  and  the 


TABLE  III.  Calculated  x-ray  structure  factors  and  directional 
compton  profiles  for  atomic  argon.  a 


AE 

Eq.  (12)  (%) 

Eq.  (5)  (?) 

S(K) 

0.2 

17,  83 

0.00 

0.05 

0.4 

17,34 

0.  00 

0. 19 

1 

14.63 

0.00 

1.02 

2 

10.19 

0.00 

2.13 

4 

6.  86 

0.00 

- 1.  69 

J(q) 

0 

5.  04 

0.00 

1.77 

0.4 

4.62 

0.00 

2.01 

1 

2.66 

0.00 

3.90 

2 

1.09 

0.00 

0.58 

4 

0.52 

0.00 

-7.08 

<r(r) 

0 

2826. 

0.02 

-0.78 

0.08 

246.  6 

-0.03 

- 1.40 

0.4 

8.61 

-1.20 

-5.38 

0.6 

1.81 

-1.60 

2.67 

1. 

0.447 

0.04 

-1.37 

1.5 

0. 186 

0.  12 

-7.27 

2 

0.0601 

0.  10 

-7.8 

4 

0.  00057 

0.06 

-7.36 

*tice  the  caption  to  Table  II  for  details. 


TABLE  IV.  Calculated  x-ray  structure  factors.  SVO,  direc- 
tional Compton  profiles,  (both  per  atom),  and  the  charge 

density  a(r)  arc  given  for  Si2.  * 


AE 

Eq.  (12)  (%) 

Eq.  (5)  (?) 

S(K)100 

0.2 

12.52 

0.01 

0.05 

0.4 

8.54 

0.03 

0.  17 

1 

-6.51 

0.09 

0.78 

2 

-2.28 

0.25 

0.78 

4 

-4.49 

0.38 

-1.35 

SW  oo, 

0.2 

13.76 

0.  00 

0.06 

0.4 

13.11 

0.01 

0.21 

1 

10.51 

0.03 

0.  80 

2 

8.13 

0.11 

0.27 

4 

5.29 

0.37 

-1.63 

«JlOo/  a^ 

0 

6.42 

0.09 

2. 16 

0.4 

3.61 

-0.33 

3.04 

1 

1.  52 

-0. 15 

-0.84 

2 

1.04 

-0. 10 

-3.  80 

4 

0.  368 

0.  65 

-1.20 

J 001/ a^ 

0 

4.79 

-0. 14 

2.07 

0.4 

3.89 

-0.  10 

2.  35 

1 

1.71 

-0.16 

0.  90 

2 

1.04 

-0.09 

-3.40 

4 

0.37 

0.  59 

-1.80 

o(r) 

0 

0.0314 

-0.24 

-6.8 

0.4435 

0.0366 

0.01 

-5.2 

1. 109 

0.129 

-0.91 

3.  09 

1.44 

0.812 

-1.42 

0.  11 

1.88 

8.39 

0.  95 

-2.37 

2.22 

1288 

0.01 

-0.44 

2.62 

6.42 

0. 12 

-2.35 

3.22 

0.218 

-1.22 

2.46 

4.22 

0.0213 

0.01 

-6.  76 

“The  column  headings  are  defined  in  the  caption  of  Table  II. 
The  line  joining  the  silicon  atoms  is  in  the  (100)  direction. 


more  correct  expression,  Eq.  (12),  are  shown  in  Tables 
n,  III,  and  IV.  In  all  cases,  the  ECP-VE  x-ray  scat- 
tering factors  and  Compton  profiles  derived  from  Eq. 
(12)  are  within  0.  7%  of  the  AE  results,  while  the  use  of 
the  approximate  density  expression,  Eq.  (5),  leads  to 
results  which  are  off  by  as  much  as  7%  for  argon  atom 
and  3%  for  the  silicon  atom  and  Si2  molecule.  The 
charge  densities  are  also  off  by  7%  for  the  approximate 
density  expression,  Eq.  (5)  and  by  less  than  1%  for  the 
more  correct  expression,  Eq.  (12). 

VI.  CONCLUSIONS 

Quantities  derivable  from  the  first-order  density  ma- 
trix can  be  reliably  obtained  from  effective-core-poten- 
tial valence-electron  calculations  to  within  roughly  one 
percent  of  the  all-electron  results  if  the  nonorthogonality 
of  the  core  and  pseudovalence  orbitals  is  correctly  taken 
into  account  (Eq.  12).  If  the  nonorthogonality  is  ignored 
and  use  is  made  of  the  standard  expression,  Eq.  (5),  er- 
rors as  large  as  7%  arc  introduced.  If  the  coreorthogo- 
nalized  valence  orbitals  arc  not  mutually  orthogonalized, 
errors  of  several  percent  are  introduced.'’ 
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APPENDIX  C 

UllF  CALCULATIONS  FOR  15-ATOM  LITHIUM  CLUSTERS 

A series  of  symmetry-restricted,  spin-unrestricted  (UHF)  Hartree- 
Fock  calculations  were  carried  out  on  a 15-atom  lithium  cluster  in  a body- 
centered  cubic  structure  with  a lattice  constant  of  3.48  angstrom  (6.57638 
bohr) . Most  of  the  calculations  were  carried  out  with  a minimal  STO-3G 
basis  set, but  a limited  study  was  also  conducted  with  a double-zeta 
basis  (obtained  by  removing  the  outermost  primitive  Gaussian  function  from 
the  2s  and  2p  contractions  in  the  ST0-3G  basis)  on  the  innermost  nine  atoms, 
while  retaining  the  fully-contracted  ST0-3G  basis  for  the  outer  six  atoms. 

As  is  seen  from  Table  I,  this  procedure  is  quite  risky,  as  it  tends  to  con- 
centrate the  charge  in  the  region  for  which  a richer  basis  set  is  provided, 
and  should  not  be  used  with  basis  sets  as  small  as  those  employed  here 

| 

(except  perhaps  for  larger  clusters,  where  the  peripheral  atoms  may  have  a 
weaker  effect  on  the  results  at  the  center  of  the  cluster). 

A multiplicity  study,  in  which  UHF  results  were  obtained  for  resultai. 
spin  S between  1/2  and  15/2  (multiplicity  2S+1  between  2 and  16)  is  described 
in  Tables  Cl  to  C4.  The  lowest  energy  is  obtained  for  S=3/2  (quartet),  un- 
like the  situation  in  diamond-structure  carbon  clusters,  which  tend  to  have 
all  peripheral  dangling-bond  electrons  with  the  same  spin.  This  appears  to 
reflect  the  metallic  structure  of  lithium,  compared  to  the  covalent  structure 
in  diamond.  However,  as  in  the  case  of  the  five-atom  carbon  cluster  (Appen- 
dices A and  D) , the  lowest  energy  is  obtained  when  there  are  no  partially- 
occupied  sets  of  degenerate  ( t or  e symmetry)  orbitals.  These  are  the  cases 
in  which  the  spin-up  and  spin-down  electrons,  separately,  form  closed-shell- 


C2 


like  structures,  while  the  higher  energy  cases  generally  involve  some  t or 
e orbitals  which  are  partially  filled  in  either  the  spin-up  or  the  spin-down 
states,  or  in  both.  On  the  other  hand,  the  most  even  charge  distributions 
are  obtained  for  S=5/2  to  9/2. 

More  useful  results,  which  could  be  directly  comparable  to  experi- 
ment, would  require  a larger  basis  set  (a  full  double-zeta  basis  on  each  atom), 
and  thus  more  computer  time  than  was  available  in  the  present  study. 


C3 


REFERENCES 


1.  W.  J.  Hehre,  R.  F.  Stewart,  and  J.  A.  Pople,  J.  Chem.  Phys.  51, 

2657  (1969). 


TABLE  Cl 


Effect 

on 


of  basis  set  in  UHF  calculations 
15-atom  lithium  cluster. 


Basis  1 

Basis  2 

Total  energy  (hartrees) 

-110.03720 

-110. 12425 

Fermi  Level  (hartrees) 

Spin  up 

-0. 1138 

-0. 1309 

Spin  down 

-0. 1040 

-0. 1092 

Charge  per  atom  (electrons) 

Shell  0 

2.8124 

3.0466 

Shell  1 

3. 0030 

3.2215 

Shell  2 

3.0273 

2.  6969 

(a)  Spin  state  S=4.  Basis  1 is  [2slp]  (ST0-3G)  on  all  15  atoms. 
Basis  2 is  [3s2p]  on  the  nine  atoms  in  shells  0 and  1,  and 
[2  sip]  on  the  six  atoms  in  shell  2. 
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FIGURE  CAPTIONS 


Fig.  Cl  Total  energy  as  a function  of  spin  state  S for  fifteen-atom 

lithium  cluster,  using  the  ST0-3G  basis. 

Fig.  C2  Highest  occupied  orbital  energy  (Fermi  level)  as  a function 

of  spin  state  S for  fifteen-atom  lithium  cluster,  using  the 
ST0-3G  basis. 

Fig.  C3  Energy  per  atom  (Mulliken  projection)  as  a function  of  spin 

state  S for  fifteen-atom  lithium  cluster,  using  ST0-3G  basis. 

Electronic  charge  per  atom  (Mulliken  gross  atomic  population) 
as  a function  of  spin  state  S for  fifteen  atom  lithium  cluster, 
using  ST0-3G  basis. 


Fig.  C4 
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USE  OF  PSEUDOPOTENTIALS  IN  CARBON 
ATOM  CLUSTER  CALCULATIONS 
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APPENDIX  D 


USE  OF  PSEUDOPOTENTIALS  IN  CARBON  ATOM  CLUSTER  CALCULATIONS 


i 


Four  series  of  symmetry-restricted,  spin-unrestricted  Hartree-Fock 
(UHF)  calculations  were  carried  out  for  the  five-atom  carbon  cluster: 

(1)  All-electron  calculations  with  a standard  double-zeta  Huzinaga- 
Dunning  basis  set.^^  These  will  be  referred  to  as  AE(SB)  (all  electron, 
standard  basis). 

(2)  All-pseudopotential  calculations,  in  which  the  two  core  elec- 

(2) 

trons  on  each  atom  were  replaced  by  a Kahn-Baybutt-Truhlar  pseudopotential, 
using  the  same  standard  basis  set  as  in  (1).  These  will  be  referred  to  as 
AP(SB)  (all  pseudo,  standard  basis). 

(3)  Pseudo-shell  calculations,  in  which  the  core  electrons  on  the 
peripheral  four  atoms  were  replaced  by  pseudopotentials,  while  all  six  electrons 
on  the  central  atom  were  treated  explicitly.  The  standard  basis  was  used. 

These  will  be  referred  to  as  PS(SB)  (pseudo  shell,  standard  basis). 

(4)  Same  as  (3),  but  the  special  valence-only  basis  of  Kahn,  Bay- 

(2) 

butt,  and  Truhlar  was  used  for  the  peripheral  four  atoms.  These  will  be 
referred  to  as  PS(VB)  (pseudo  shell,  valence  basis). 

A comparison  of  series  (1)  and  (2)  allows  a general  assessment  of  the 
effect  of  the  use  of  pseudopotentials  on  molecular  computational  results.  A 
comparison  of  series  (2)  and  (3)  allows  an  examination  of  the  idea  of  treating 
peripheral  atoms  in  a cluster  by  a lower-order  approximation  than  was  used  for 
the  central  atom  (or  atoms).  This  examination  is  facilitated  by  examining 

(3  4) 

properties  projected  for  the  central  atom  in  the  cluster.  ' A comparison  of 
series  (3)  and  (4)  provides  an  assessment  of  the  effect  of  the  use  of  a valence- 
only  basis  for  those  atoms  which  are  treated  in  a pseudopotential  approximation. 
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The  use  of  such  a basis  is  essential  if  the  full  economies  of  the  pseudo- 
potential approach  are  to  be  realized  in  calculations  on  larger  clusters. 

For  each  series,  calculations  were  carried  out  for  total  spin  S 

ranging  from  zero  to  six  or  seven  (multiplicity  one  to  thirteen  or  fifteen), 

(3) 

and  as  in  previous  calculations,  it  was  found  that  the  results  for  S=6, 
with  twelve  parallel  spins  representing  the  twelve  dangling  bonds  on  the 
peripheral  four  atoms,  produced  the  lowest  cluster  energy  and  most  closely 
represented  the  situation  in  the  diamond  cluster. 

The  principal  results  of  these  calculations  are  shown  in  Tables  DI 
to  DV  and  Figures  Dl  to  D4.  A more  complete  account  is  being  prepared  for 
publication. 

It  is  clear  from  the  results  given  that  the  calculations  which 
utilize  the  core  approximation  reproduce  the  trends  and  characteristics  of 
the  all-electron  calculations.  Central-atom  (projected)  results  for  the  pseudo- 
shell (PS)  calculations  are  particularly  close  to  the  all-electron  results, 
while  totai  cluster  results  for  all  the  pseudopotential  calculations  are 
quite  close  to  each  other,  and  often  could  not  be  resolved  from  each  other 
on  the  scale  of  the  plots.  For  the  variation  of  energy  with  lattice  constant, 
the  pseudo-shell  results  are  seen  to  reproduce  the  all-electron  results  much 
better  than  the  all-pseudo  treatments. 

The  use  of  a valence  basis  on  those  atoms  for  which  the  core  approxi- 
mation is  invoked  is  seen  to  have  little  effect  on  the  results  except  in  the 
case  of  the  vacant  orbital  energies.  As  is  well  known  from  molecular  calcu- 
lations, the  vacant  orbitals  are  particularly  sensitive  to  the  choice  of  basis, 
and  have  little  individual  significance  except  in  minimal  basis  calculations. 
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Thus  it  is  seen  that  the  use  of  a core  approximation  together  with  a valence 
basis  is  a viable  approach  for  the  treatment  of  peripheral  atoms  in  cluster 
calculations.  In  treatments  of  clusters  of  heavier  atoms,  successively  more 
electron  shells  can  be  treated  by  this  approximation  as  one  proceeds  from  the 
center  of  the  cluster  to  its  periphery,  resulting  in  significant  economies  in 


computational  effort. 
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TABLE  DI. 

Energies  of  five-atom  carbon  cluster  as  a function  of  spin  state  S, 
for  different  types  of  calculation,  in  hartrees. 


s 

Total 

Energy 

Central 

Atom 

Energy 

AE(SB) 

AE(SB) 

AP (SB) 

PS (SB) 

PS(VB) 

0 

-188.4019 

-188.4632 

-188.4549 

-188.4580 

-37.9080 

1 

.3030 

.3643 

.3562 

.3602 

.9042 

2 

.4520 

.5142 

.5059 

.5114 

.9144 

3 

.3690 

.4293 

.4209 

.4289 

.8927 

4 

.2452 

.3245 

.3157 

.3285 

.8363 

5 

.4076 

.4758 

.4664 

.4844 

.8852 

6 

.6380 

.7122 

.7055 

.7183 

.9291 

7 

-187.8963 

-187.9959 

-187.9733 

- 

.6102 

(a) 

The  notation  for  the  different  types  of  calculation  is  described  in 
the  text.  The  energies  for  the  pseudopotential  calculations  include 
core  contributions  computed  from  the  difference  in  the  energies  of  single- 
atom calculations  with  and  without  the  core  approximation  (with  the 
corresponding  bases).  The  lattice  constant  is  3.56  angstrom. 
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TABLE  DII. 

Ionization  potentials  of  five-atom  cluster  as  a 

function  of  spin  state  S,  for  different  types  of 

(a) 

calculation,  in  electron  volts. 


s 

AE (SB) 

AP (SB) 

PS  (SB) 

PS(VB) 

0 

7.96 

4 

8.11 

4 

8.08 

+ 

8.18  + 

1 

4.65 

4 

4.05 

4 

4.04 

4.14  1 

2 

7.31 

4 

7.45 

4 

7.45 

7.57  + 

3 

‘ 9.07 

4 

9.16 

4 

9.16 

9.10  4 

4 

5.35 

4 

5.05 

4 

5.05 

+ 

4.98  4 

5 

7.50 

4 

7.47 

4 

7.36 

7.25  4 

6 

11.01 

4 

11.10 

+ 

11.07 

+ 

11.11  4 

7 

-0.02 

4 

0.56 

0.02 

t 

- 

Csl') 

The  ionization  potentials  are  the  Koopmans 
theorem  values  (i.e.,  the  negatives  of  the 
highest  occupied  orbital  energy).  The  arrows 
next  to  the  entries  indicate  whether  the  high- 
est occupied  orbital  has  spin  up  or  spin  down. 


TABLE  Dill. 


dompton  profiles  (per  atom)  of  five-atom  cluster 
in  (1,0,0)  direction  for  S=6,  for  different  types 

(3) 

of  calculation. 


q 

AE (SB) 

AP (SB) 

PS (SB) 

PS(VB) 

For  total 

cluster 

0.0 

2.464 

2.508 

2.498 

2.458 

0.5 

1.985 

1.975 

1.976 

1.989 

1.0 

1.049 

1.029 

1.034 

0.981 

1.5 

0.504 

0.496 

0.498 

0.500 

2.0 

0.303 

0.302 

0.303 

0.304 

For  atom  1 

0.0 

2.354 

2.446 

2.382 

2.517 

0.5 

1.964 

1.972 

1.977 

2.097 

1.0 

1.049 

1.026 

1.050 

1.071 

1.5 

0.486 

0.474 

0.484 

0.486 

2.0 

0.289 

0.287 

0.288 

0 . 288 

For  atoms 

2-5 

0.0 

2.492 

2.524 

2.526 

2.444 

0.5 

1.990 

1.975 

1.976 

1.962 

1.0 

1.050 

1.030 

1.030 

1.030 

1.5 

0.508 

0.501 

0.501 

0.504 

2.0 

0.306 

0.306 

0.306 

0.307 

Atom  1 is  the  central  atom,  in  the  (0,0,0) 
position.  The  core  contribution  to  the  pseudo- 
potential results  was  taken  from  a comparison  of 
single-atom  calculations  with  and  without  the 
core  approximation. 
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TABLE  DIV 

Compton  profiles  (per  atom)  of  five-atom  cluster 

in  (1,1,1)  direction  for  S=6,  for  different  types 
C 3 ) 

of  calculation. 


q 

AE(SB) 

AP  (SB) 

PS (SB) 

PS(VB) 

For  total 

cluster 

0.0 

2.431 

2.476 

2.465 

2.427 

0.5 

1.976 

1.968 

1.969 

1.979 

1.0 

1.080 

1.062 

1.066 

1.070 

1.5 

0.508 

0.514 

0.499 

0.503 

2.0 

0.295 

0.293 

0.294 

0.295 

For  atom 

1 

0.0 

2.315 

2.405 

2.341 

2.480 

0.5 

1.927 

1.941 

1.944 

2.062 

1.0 

1.120 

1.100 

1.124 

1.144 

1.5 

0.496 

0.474 

0.489 

0.492 

2.0 

0.271 

0.264 

0.267 

0.267 

For  atom 

2 

0.0 

2.573 

2.601 

2.602 

2.558 

0.5 

1.942 

1.937 

1.934 

1.906 

1.0 

1.128 

1.113 

1.112 

1.102 

1.5 

0.491 

0.472 

0.477 

0.482 

2.0 

0.269 

0.263 

0.265 

0.268 

For  atoms 

3-5 

0.0 

2.423 

2.457 

2.460 

2.366 

0.5 

2.004 

1.987 

1.989 

1.976 

1.0 

1.051 

1.032 

1.031 

1.034 

1.5 

0.517 

0.511 

0.510 

0.513 

2.0 

0.312 

0.312 

0.312 

0.313 

Atom  1 is  the  central  atom  and  atom  2 is  in 
the  (1,1,1)  position.  The  core  contribution  to 
the  pseudopotential  result  was  taken  from  a 
comparison  of  single  atom  calculations  with  and 
without  the  core  approximation. 
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TABLE  DV. 

Compton  profile  of  the  central  atom  in  five  atom 

cluster  at  q=0  as  a function  of  spin  state  S, 

(a) 

for  different  types  of  calculation. 


s 

AE(SB) 

AP (SB) 

PS (SB) 

PS(VB) 

(1,0,0) 

direction 

0 

2.117 

2.188 

2.115 

2.145 

1 

2.181 

2.258 

2.184 

2.182 

2 

2.281 

2.357 

2.286 

2.269 

3 

2.302 

2.382 

2.309 

2.295 

4 

2.187 

2.288 

2.215 

2.213 

5 

2.124 

2.209 

2.137 

2.190 

6 

2.354 

2.446 

2.382 

2.517 

7 

1.883 

2.075 

1.960 

- 

(1,1,1) 

direction 

0 

2.341 

2.424 

2.354 

2.401 

1 

2.386 

2.470 

2.399 

2.411 

2 

2.470 

2.554 

2.484 

2.479 

3 

2.428 

2.507 

2.435 

2.436 

4 

2.223 

2.307 

2.236 

2.253 

5 

2.041 

2.120 

2.050 

2.119 

6 

2.315 

2.405 

2.341 

2.480 

7 

2.436 

2.570 

2.504 

(a) 

' 'The  core  contribution  to  the  AP(SB)  result 
was  obtained  from  a comparison  of  single-atom 
calculations  with  and  without  the  core  approxi- 
mation. 


DIO 


FIGURE  CAPTIONS 

Fig.  D1  Total  energy  of  five-atom  cluster  as  function  of  spin  state 

S,  for  different  types  of  calculation.  The  notation  is 
described  in  the  text.  The  results  for  PS(VB)  are  indistin- 
guishable from  those  for  PS(SB)  on  the  scale  of  this  graph. 

See  also  Table  I. 

Fig.  D2  Orbital  energies  of  five-atom  cluster  for  S=6,  for  different 

types  of  calculations.  The  numbers  next  to  the  orbital  levels 
specify  the  degeneracy  of  each  orbital. 

Fig.  D3  Valence  band  width  and  band  gap  of  five-atom  cluster  as  a 

function  of  spin  state  S,  for  different  types  of  calculation. 
The  results  for  PS(SB)  and  PS(VB)  are  indistinguishable  from 
those  for  SP(SB)  on  the  scale  of  this  figure.  The  points 
marked  with  an  asterisk  are  those  for  which  there  is  overlap 
between  an  occupied  spin-up  band  and  a vacant  spin-down  band, 
or  vice  versa.  This  happens  when  a set  of  degenerate  orbitals 
is  partially  filled,  and  results  from  the  imposition  of  spatial 
symmetry  on  the  orbitals. 

Fig.  DA  Cluster  energy  of  five-atom  cluster  as  a function  of  lattice 

constant,  for  S=6,  for  different  types  of  calculation. 


Fig.  D3 
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